      subroutine sseig

      use param
      use allocs

#include <global.h>
      include 'ntbytes.h'
      include 'geom.h'
      include 'geomf.h'
      include 'xsec.h'
      include 'sanm.h'
      include 'sanmjin.h'
      include 'nem.h'
      include 'setls.h'
      include 'eigv.h'
      include 'cmfd.h'
      include 'cmfdp1.h'
      include 'scmfdp1.h'

      character*30 t1
      real(NBF) const1,const2

! initial condition
      eigv=1
      reigv=1/eigv
#define init
#ifdef init
      do im=1,ng
        do i=1,tnode
          avgflx(1,i,im)=1 
          avgflx(2,i,im)=0.01
          psphi(2,i,im)=avgflx(2,i,im)
          psphi(1,i,im)=avgflx(1,i,im)+2*psphi(2,i,im)
          qs(0,i,im)=1 
!          sfrt(i,im)=1
!          jinl(1,i,im)=0.6
!          jinr(1,i,im)=0.5
!          jinl(2,i,im)=-0.01
!          jinr(2,i,im)=-0.05
          ajil(1,i,im)=0.5*jnet(1,i,im)+0.25*avgflx(1,i,im)+0.3125*avgflx(2,i,im)
          ajir(1,i,im)=-0.5*jnet(1,i+1,im)+0.25*avgflx(1,i,im)+0.3125*avgflx(2,i,im)
          ajil(2,i,im)=0.5*jnet(2,i,im)-0.0625*avgflx(1,i,im)+0.3125*avgflx(2,i,im)
          ajir(2,i,im)=-0.5*jnet(2,i+1,im)-0.0625*avgflx(1,i,im)+0.3125*avgflx(2,i,im)
        enddo
!
        do j=1,fnode
          aflx(1,j,im)=1 
          aflx(2,j,im)=0.01
          pphi(2,j,im)=aflx(2,j,im)
          pphi(1,j,im)=aflx(1,j,im)+2*pphi(2,j,im)
          qsf(0,j,im)=1 
!          jinl(1,j,im)=0.5*jnetf(1,j,im)+0.25*aflx(1,j,im)+0.3125*aflx(2,j,im)
!          jinr(1,j,im)=-0.5*jnetf(1,j+1,im)+0.25*aflx(1,j,im)+0.3125*aflx(2,j,im)
!          jinl(2,j,im)=0.5*jnetf(2,j,im)-0.0625*aflx(1,j,im)+0.3125*aflx(2,j,im)
!          jinr(2,j,im)=-0.5*jnetf(2,j+1,im)-0.0625*aflx(1,j,im)+0.3125*aflx(2,j,im)        
        enddo
        call dumflx
      enddo

      do j=1,tnode
        sum=0
        do im=1,ng
          sum=sum+xsnf(j,im)*avgflx(1,j,im)
        enddo
        psif(j)=sum*hz(j)
      enddo
#endif

! Ready
! transfer group/composition xsec to region/node xsec
      do im=1,ng
        do i=1,tnode
          sigd(1,i,im)=4*xsd(1,i,im)*(rhz(i)**2) ! SIGD0
          sigd(2,i,im)=4*xsd(2,i,im)*(rhz(i)**2) ! SIGD2 
          rsigd(:,i,im)=1/sigd(:,i,im)
        enddo
!
        do j=1,fnode
          sigdf(1,j,im)=4*xsdf(1,j,im)*(rhzf(j)**2) ! SIGD0
          sigdf(2,j,im)=4*xsdf(2,j,im)*(rhzf(j)**2) ! SIGD2 
          rsigdf(:,j,im)=1/sigdf(:,j,im)
        enddo
      enddo

! constants for the SANM    
      do im=1,ng
        do i=1,tnode
          const1=sqrt((3*sigd(2,i,im)*xsr(i,im))**2+6*sigd(1,i,im)*sigd(2,i,im)*xsr(i,im) &
                *(4*xsr(i,im)-5*xst(i,im))+(sigd(1,i,im)*(4*xsr(i,im)+5*xst(i,im)))**2)
          const2=sqrt(-60*sigd(1,i,im)*sigd(2,i,im)*xsr(i,im)*xst(i,im) &
                +(3*sigd(2,i,im)*xsr(i,im)+sigd(1,i,im)*(4*xsr(i,im)+5*xst(i,im)))**2)
! eigenvectors
          eigvec1(i,im)=-(4*sigd(1,i,im)*xsr(i,im)+3*sigd(2,i,im)*xsr(i,im)-5*sigd(1,i,im)*xst(i,im) &
                       -const2)/(4*sigd(1,i,im)*xsr(i,im))
          eigvec2(i,im)=-(4*sigd(1,i,im)*xsr(i,im)+3*sigd(2,i,im)*xsr(i,im)-5*sigd(1,i,im)*xst(i,im) &
                       +const2)/(4*sigd(1,i,im)*xsr(i,im))
! isiD
          imq1(i,im)=(sigd(1,i,im)*(4*xsr(i,im)+5*xst(i,im))+3*sigd(2,i,im)*xsr(i,im)-const2) &
                    /(3*sigd(2,i,im)*const2)
          imq2(i,im)=(-sigd(1,i,im)*(4*xsr(i,im)+5*xst(i,im))-3*sigd(2,i,im)*xsr(i,im)-const2) &
                    /(3*sigd(2,i,im)*const2)
! eigenvalues
          ksq(1,i,im)=(-(sigd(1,i,im)*(4*xsr(i,im)+5*xst(i,im)))**2+3*sigd(2,i,im) &
                     *xsr(i,im)*(-3*sigd(2,i,im)*xsr(i,im)+const1)+sigd(1,i,im) &
                     *(-6*sigd(2,i,im)*xsr(i,im)*(4*xsr(i,im)-5*xst(i,im))+(4*xsr(i,im)+5*xst(i,im))*const1)) &
                     /(6*sigd(1,i,im)*sigd(2,i,im)*const2)
         ksq(2,i,im)=((sigd(1,i,im)*(4*xsr(i,im)+5*xst(i,im)))**2+3*sigd(2,i,im) &
                     *xsr(i,im)*(3*sigd(2,i,im)*xsr(i,im)+const1)+sigd(1,i,im) &
                     *(6*sigd(2,i,im)*xsr(i,im)*(4*xsr(i,im)-5*xst(i,im))+(4*xsr(i,im)+5*xst(i,im))*const1)) &
                     /(6*sigd(1,i,im)*sigd(2,i,im)*const2) 
          krt(:,i,im)=sqrt(ksq(:,i,im))
          kqu(:,i,im)=ksq(:,i,im)**2
          ksx(:,i,im)=ksq(:,i,im)**3
          rkrt(:,i,im)=1/krt(:,i,im)
          rksq(:,i,im)=1/ksq(:,i,im)
          rkqu(:,i,im)=1/kqu(:,i,im)
          rksx(:,i,im)=1/ksx(:,i,im)
        enddo
      enddo

! surface flux 
      do im=1,ng
        do j=1,tnode
          jm1=j-1
          jp1=j+1
          btm1=beta(1,jm1,im)
          bt1=beta(1,j,im)
          btp1=beta(1,jp1,im)
          btm2=beta(2,jm1,im)
          bt2=beta(2,j,im)
          sflux(1,j,im)=(bt1*avgflx(1,j,im)+btm1*avgflx(1,jm1,im))/(bt1+btm1)  ! surface flux for the FDM
          sflux(2,j,im)=(bt2*avgflx(2,j,im)+btm2*avgflx(2,jm1,im))/(bt2+btm2)
        enddo
        sflux(1,j,im)=avgflx(1,j-1,im)
        sflux(2,j,im)=avgflx(2,j-1,im)
      enddo
! Start SANM SP3
!      t1="dprob.txt"
!      open (unit=84,file=t1,status="replace")
      do iout=1,nout
        if (opt.eq.0) then
          call cmfd
          call solsanm
        elseif (opt.eq.1) then
!          call cmfdp1
          call scmfdp1
          call solsanm
          if (iout.lt.50) then
            do im=1,ng
              do j=1,tnode+1
                sfrt(j,im)=0
              enddo
            enddo
          endif
        elseif (opt.eq.2) then
!          call cmfdp1
          call scmfdp1
!          call cmfd
          call trangl
          call solnem
          call tranlg   
          continue   
        endif
! check a dhat problem
!        write(84,*) iout,2,2,dprob(2,2) 

        err2d=1
        err2=0  
        econv=0
        psipsi=0
        psipsid=0

        do j=1,tnode
          psifd(j)=psif(j)
        enddo

        do j=1,tnode
          sum=0
          do im=1,ng
            sum=sum+xsnf(j,im)*avgflx(1,j,im)
          enddo
          psif(j)=sum*hz(j)
          psipsi=psipsi+psif(j)*psif(j)
          psipsid=psipsid+psif(j)*psifd(j)
          err=psifd(j)-psif(j)
          err2=err2+err*err
        enddo

        domr=sqrt(err2/err2d)
!
        err2d=err2
        eigvd=eigv
        eigv=eigv*psipsi/psipsid5
        reigv=1/eigv  

        print*,iout,eigv,abs(eigv-eigvd)/abs(eigv)

        econv=sqrt(err2/psipsid)

!     exit condition
        if (abs(eigv-eigvd).lt.kconv .and. econv.lt.fconv) then 
          niter=iout
          print*,maxval(pn),"box"
          exit
        endif
      enddo

!#define test1
#ifdef test1

      allocate (tavg1(nreg,ng))
      allocate (tavg2(nreg,ng))

      do im=1,ng
        ib=1
        ie=0
        do j=1,nreg
          ie=ie+sn(j)
          temp1=0
          temp2=0
          do i=ib,ie
            temp1=temp1+avgflx(1,i,im)
            temp2=temp2+avgflx(2,i,im)
          enddo
          tavg1(j,im)=temp1/sn(j)
          tavg2(j,im)=temp2/sn(j)
          ib=i
        enddo
      enddo
#endif

      end subroutine
